###############################################################################
############# Conjoint Experiment BJPolS  ##############
###############################################################################

# load in packages
library(dplyr)
library(tidyr)
library(forcats)
library(cregg)
library(foreign)
library(data.table) 
library(ggplot2)
library(psych) 
library(sjstats)
library(sjmisc)
library(haven)
library(tidyverse)
library(broom)
library(hrbrthemes)
library(egg)
library(grid)


# load in data
cj_data <- read.csv("conjointdata_BJPolS.csv")


######### Preparing data for conjoint #########################################

#### cj factorial variables ####
cj_data$issue<-factor(cj_data$issue, labels = c("Emissions", "Refugees" , "Currency" , "Foreign aid"))
cj_data$initiative <-factor(cj_data$initiative, labels = c("NGO", "Government"))
cj_data$recruitment <-factor(cj_data$recruitment, labels = c("Random selection", "Self-selection"))
cj_data$size <-factor(cj_data$size, labels = c("Small", "Medium" , "Large"))
cj_data$composition <-factor(cj_data$composition, labels = c("Citizens alone", "Mixed groups"))
cj_data$output <-factor(cj_data$output, labels = c("In favor of measure", "Against measure"))
cj_data$consensus <-factor(cj_data$consensus, labels = c("Narrow majority", "Clear majority"))
cj_data$format <-factor(cj_data$format, labels = c("Face-to-face", "Online"))
cj_data$authorization <-factor(cj_data$authorization, labels = c("Recommendation", "Referendum" , "Binding decision"))
cj_data$outcome.favorability <- factor(cj_data$outcome.fav, labels = c("Preference mismatch", "Preference match"))

feature_names <- c(
  'issue' = "Policy Issue",
  'initiative' = "Initiative",
  'recruitment' = "Recruitment",
  'size' = "Group Size",
  'composition' = "Group Composition",
  'consensus' = "Degree of Consensus",
  'format' = "Discussion Format",
  'authorization' = "Authorization",
  'outcome.favorability' = "Outcome Favorability"
)

#### create timing variable for conjoint screens (> median engagement with scenarios) ####
describe(cj_data$page_p_screen1_timing)
describe(cj_data$page_p_screen2_timing)
describe(cj_data$page_p_screen3_timing)
describe(cj_data$page_p_screen4_timing)
describe(cj_data$page_p_screen5_timing)
describe(cj_data$page_p_screen6_timing)
cj_data$select_timing1 <- ifelse(cj_data$comparision_table == 1 & cj_data$page_p_screen1_timing > 49.98, 1, 0)
cj_data$select_timing2 <- ifelse(cj_data$comparision_table == 2 & cj_data$page_p_screen2_timing > 32.15, 1, 0)
cj_data$select_timing3 <- ifelse(cj_data$comparision_table == 3 & cj_data$page_p_screen3_timing > 28.42, 1, 0)
cj_data$select_timing4 <- ifelse(cj_data$comparision_table == 4 & cj_data$page_p_screen4_timing > 26.01, 1, 0)
cj_data$select_timing5 <- ifelse(cj_data$comparision_table == 5 & cj_data$page_p_screen5_timing > 23.97, 1, 0)
cj_data$select_timing6 <- ifelse(cj_data$comparision_table == 6 & cj_data$page_p_screen6_timing > 22.74, 1, 0)
cj_data$timing <- rowSums(subset(cj_data, select = c(select_timing1, select_timing2, select_timing3, select_timing4, select_timing5, select_timing6)), na.rm = TRUE)
cj_data$timing <- factor(cj_data$timing, labels = c("below median", "above median"))


############################################################################### 
######### MAIN ################################################################

######### Estimating AMCEs using Cregg package (Leeper 2020) ##############
#### specifying model ####
Mchoice <- choice ~ issue + initiative + recruitment + size + composition + consensus + format + authorization + outcome.favorability # conjoint with all attributes, outcome favorability, choice outcome

#### baseline models ####
amce_choice <- cj(cj_data, Mchoice, id = ~ID,  estimate="amce", feature_labels = feature_names)

#### Plot AMCE choice ####
plot(amce_choice, feature_headers = FALSE)+
  ggplot2::facet_grid(feature~.,  scales = "free_y", space = "free", switch ="y") +
  theme(panel.border = element_rect(linetype = "blank"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 12),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none",
        panel.spacing.y = unit(1, "lines"),
        strip.background= element_rect(color="white", size = 2, fill = "#f9f9f9"), 
        strip.text.y.left = element_text(angle = 0),
        strip.placement = "outside",
        strip.text.y = element_text(colour = "#091926", size = 12, face = "bold")) +
  geom_point(size=2) +
  scale_color_manual(values = rep("#062f51", 9))

######### Disaffected citizens ################################################
cj_satisfaction <- cj_data %>% select(ID:authorization, q2:q3_3, q6_1, q6_2, outcome.favorability) # subset including all relevant variables

#### 1. Political dissatisfaction ####
# q2 (satisfaction democracy)

sub_satisfaction <- cj_data %>% select(ID:authorization, q2, outcome.favorability) %>% drop_na()
f1 <- table(sub_satisfaction$q2)
f_sat <- prop.table(f1)
f_sat <-round(f_sat, digits = 3)
f1kum <- cumsum(f1)
f_satkum <- cumsum(f_sat)
f_satkum <- round(f_satkum, digits = 3)
cbind(f1, f_sat, f1kum, f_satkum)

# median split
describe(sub_satisfaction$q2)
sub_satisfaction$satsplit[sub_satisfaction$q2 >= 5] <- 1L
sub_satisfaction$satsplit[sub_satisfaction$q2 < 5] <- 2L
sub_satisfaction$satsplit <- factor(sub_satisfaction$satsplit, 1:2, c("High satisfaction", "Low satisfaction"))
describeBy(evaluation ~ satsplit,data =sub_satisfaction)
# differences in means between satisfied vs dissatisfied citizens?
t.test(evaluation~satsplit, var.equal = TRUE, alternative = "two.sided", data = sub_satisfaction) # t=9.95, p<0.001
# nested model comparison test (does any of the interactions between the "by variable" and feature levels differ from zero?)
cj_anova(sub_satisfaction, Mchoice, by = ~satsplit)

#### 2. External efficacy ####
# q6_1 to q6_2 (external efficacy)

alpha(subset(cj_satisfaction, select = c(q6_1, q6_2)), check.keys = TRUE) # alpha .86
cj_satisfaction$external_index <- rowMeans(subset(cj_satisfaction, select = c(q6_1, q6_2)), na.rm = TRUE)
sub_external <- cj_satisfaction %>% filter(external_index != "NaN")
f2 <- table(cj_satisfaction$external_index)
f_ext <- prop.table(f2)
f_ext <-round(f_ext, digits = 3)
f2kum <- cumsum(f2)
f_extkum <- cumsum(f_ext)
f_extkum <- round(f_extkum, digits = 3)
cbind(f2, f_ext, f2kum, f_extkum)

# median split
describe(cj_satisfaction$external_index)
sub_external$externalsplit[sub_external$external_index >= 2.5] <- 1L
sub_external$externalsplit[sub_external$external_index < 2.5] <- 2L
sub_external$externalsplit <- factor(sub_external$externalsplit, 1:2, c("High ext. efficacy", "Low ext. efficacy"))
cj_anova(sub_external, Mchoice, by = ~externalsplit)
describeBy(evaluation ~ externalsplit,data =sub_external)
t.test(evaluation~externalsplit, var.equal = TRUE, alternative = "two.sided", data = sub_external) # t=12.3, p<0.001
cj_anova(sub_external, Mchoice, by = ~externalsplit)

#### 3. Stealth  ####
#  q13_1:q13_4 
cj_stealthsun <- cj_data %>% select(ID:authorization, q7_1:q7_4, q13_1:q13_4, outcome.favorability)

alpha(subset(cj_stealthsun, select = c(q13_1:q13_4)), check.keys = TRUE) # alpha .66
cj_stealthsun$stealth_index <- rowMeans(subset(cj_stealthsun, select = c(q13_1, q13_2, q13_3, q13_4)), na.rm = TRUE)
sub_stealth <- cj_stealthsun %>% filter(stealth_index != "NaN")
f3 <- table(sub_stealth$stealth_index)
f_stealth <- prop.table(f3)
f_stealth <-round(f_stealth, digits = 3)
f3kum <- cumsum(f3)
f_stealthkum <- cumsum(f_stealth)
f_stealthkum <- round(f_stealthkum, digits = 3)
cbind(f3, f_stealth, f3kum, f_stealthkum)

# median split
describe(sub_stealth$stealth_index)
sub_stealth$stealthsplit[sub_stealth$stealth_index <= 4] <- 1L
sub_stealth$stealthsplit[sub_stealth$stealth_index > 4] <- 2L
sub_stealth$stealthsplit <- factor(sub_stealth$stealthsplit, 1:2, c("Low stealth", "High stealth"))
describeBy(evaluation ~ stealthsplit,data =sub_stealth)
t.test(evaluation~stealthsplit, var.equal = TRUE, alternative = "two.sided", data = sub_stealth) # t=9.23, p<0.001
cj_anova(sub_stealth, Mchoice, by = ~stealthsplit)

#### 4. Populism ####
cj_pop <- cj_data %>% select(ID:authorization, q12_1:q12_7,q13_1, outcome.favorability)

alpha(subset(cj_pop, select = c(q12_1:q13_1)), check.keys = TRUE) # alpha .81
alpha(subset(cj_pop, select = c(q12_1, q12_3, q12_7)), check.keys = TRUE) # 1. dimension homogeneity # alpha .86
alpha(subset(cj_pop, select = c(q12_2, q12_5, q12_6)), check.keys = TRUE) # 2. dimension sovereignty # alpha .81
alpha(subset(cj_pop, select = c(q12_4, q13_1)), check.keys = TRUE) # 3. dimension anti elitism # alpha .61
alpha(subset(cj_pop, select = c(q12_2, q12_4, q12_5, q12_6, q13_1)), check.keys = TRUE) # 2+3 dimension # alpha .8

cj_pop$pop_index <- rowMeans(subset(cj_pop, select = c(q12_1, q12_2, q12_3, q12_4, q12_5, q12_6, q12_7,q13_1)), na.rm = TRUE)
sub_pop <- cj_pop %>% filter(pop_index != "NaN")
f4 <- table(cj_pop$pop_index)
f_pop <- prop.table(f4)
f_pop <-round(f_pop, digits = 3)
f4kum <- cumsum(f4)
f_popkum <- cumsum(f_pop)
f_popkum <- round(f_popkum, digits = 3)
cbind(f4, f_pop, f4kum, f_popkum)

# median split
describe(sub_pop$pop_index)
sub_pop$popsplit[sub_pop$pop_index <= 4.62] <- 1L
sub_pop$popsplit[sub_pop$pop_index > 4.62] <- 2L
sub_pop$popsplit <- factor(sub_pop$popsplit, 1:2, c("Low populist", "High populist"))
describeBy(evaluation ~ popsplit,data =sub_pop)
t.test(evaluation~popsplit, var.equal = TRUE, alternative = "two.sided", data = sub_pop) # t=12.93, p<0.001
cj_anova(sub_pop, Mchoice, by = ~popsplit)

#### PLOTS Subgroups ####
diff_amce_sat_choice <- cj(sub_satisfaction, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~satsplit)
diff_amce_external_choice <- cj(sub_external, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~externalsplit)
diff_amce_stealth_choice <- cj(sub_stealth, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~stealthsplit)
diff_amce_pop_choice <- cj(sub_pop, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~popsplit)

diff_disent_estimates_choice <-bind_rows(
  diff_amce_sat_choice %>% mutate(model = "Dissatisfied citizens"),
  diff_amce_external_choice %>% mutate(model = "Low-efficacious citizens"),
  diff_amce_stealth_choice %>% mutate(model = "Stealth citizens"),
  diff_amce_pop_choice %>% mutate(model = "Populist citizens"))

diff_disent_choice <- diff_disent_estimates_choice %>% select(BY:level, estimate:model)
diff_disent_choice$model <- factor(diff_disent_choice$model, levels = c("Dissatisfied citizens", "Low-efficacious citizens", "Stealth citizens", "Populist citizens"))

plot(diff_disent_choice, feature_headers = FALSE, xlab = "Estimated Difference (AMCE) for Disaffected Citizens") + 
  ggplot2::facet_wrap(~model, ncol = 4) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 10),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = rep("#833d3d", 9))


############################################################################### 
######### APPENDIX ############################################################

######### Benchmark models ####################################################
Mrating <- evaluation ~ issue + initiative + recruitment + size + composition + consensus + format + authorization + outcome.favorability # conjoint with all attributes, outcome favorability, rating outcome
Mchoice <- choice ~ issue + initiative + recruitment + size + composition + consensus + format + authorization + outcome.favorability # conjoint with all attributes, outcome favorability, rating outcome

mm_choice <- cj(cj_data, Mchoice, id = ~ID, estimate="mm", h0=0.5, feature_labels = feature_names)
mm_rating <- cj(cj_data, Mrating, id = ~ID, estimate="mm", feature_labels = feature_names)
amce_choice<- cj(cj_data, Mchoice, id = ~ID, estimate="amce", feature_labels = feature_names)
amce_rating<- cj(cj_data, Mrating, id = ~ID, estimate="amce", feature_labels = feature_names)

#### Figure B.1 benchmark model rating AMCE ####
plot(amce_rating, feature_headers = FALSE)+
  ggplot2::facet_grid(feature~.,  scales = "free_y", space = "free", switch ="y") +
  theme(panel.border = element_rect(linetype = "blank"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 12),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none",
        panel.spacing.y = unit(1, "lines"),
        strip.background= element_rect(color="white", size = 2, fill = "#f9f9f9"), 
        strip.text.y.left = element_text(angle = 0),
        strip.placement = "outside",
        strip.text.y = element_text(colour = "#091926", size = 12, face = "bold")) +
  geom_point(size=2) +
  scale_color_manual(values = rep("#062f51", 9))

######### Subgroups ###########################################################
##### 1. Combined difference plots #####
#### Figure B2.1: Combined difference plots AMCE vs MM for satisfaction ####
diff_amce_sat_rating <- cj(sub_satisfaction, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~satsplit)
diff_mm_sat_rating <- cj(sub_satisfaction, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~satsplit)
diff_amce_sat_choice <- cj(sub_satisfaction, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~satsplit)
diff_mm_sat_choice <- cj(sub_satisfaction, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names,feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~satsplit)
diff_amce_sat_rating$Estimate <- "AMCE (Rating)"
diff_mm_sat_rating$Estimate <- "Marginal Means (Rating)"
diff_amce_sat_choice$Estimate <- "AMCE (Choice)"
diff_mm_sat_choice$Estimate <- "Marginal Means (Choice)"
diff_sat_rating <- rbind(diff_amce_sat_rating, diff_mm_sat_rating)
diff_sat_choice <- rbind(diff_amce_sat_choice, diff_mm_sat_choice)

sat1 <- plot(diff_sat_rating, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

sat2 <- plot(diff_sat_choice, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

combined_sat <- ggarrange(sat1 + theme(axis.ticks.y = element_blank(),
                                       plot.margin = margin(r = 1) ),
                          sat2 + theme(axis.text.y = element_blank(),
                                       axis.ticks.y = element_blank(),
                                       axis.title.y = element_blank(),
                                       plot.margin = margin(l = 1)),
                          nrow = 1,
                          bottom = textGrob("Estimated differences for 'political dissatisfaction' compared to 'political satisfaction'", vjust = 0.5, gp = gpar(col="#091926", fontsize=12)))

#### Figure B2.2: Combined difference plots AMCE vs MM for external efficacy ####
diff_amce_external_rating <- cj(sub_external, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~externalsplit)
diff_mm_external_rating <- cj(sub_external, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~externalsplit)
diff_amce_external_choice <- cj(sub_external, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~externalsplit)
diff_mm_external_choice <- cj(sub_external, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~externalsplit)
diff_amce_external_rating$Estimate <- "AMCE (Rating)"
diff_mm_external_rating$Estimate <- "Marginal Means (Rating)"
diff_amce_external_choice$Estimate <- "AMCE (Choice)"
diff_mm_external_choice$Estimate <- "Marginal Means (Choice)"
diff_external_rating <- rbind(diff_amce_external_rating, diff_mm_external_rating)
diff_external_choice <- rbind(diff_amce_external_choice, diff_mm_external_choice)

external1 <- plot(diff_external_rating, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

external2 <- plot(diff_external_choice, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

combined_external <- ggarrange(external1 + theme(axis.ticks.y = element_blank(),
                                                 plot.margin = margin(r = 1) ),
                               external2 + theme(axis.text.y = element_blank(),
                                                 axis.ticks.y = element_blank(),
                                                 axis.title.y = element_blank(),
                                                 plot.margin = margin(l = 1)),
                               nrow = 1,
                               bottom = textGrob("Estimated differences for 'low external efficacy' compared to 'high external efficacy'", vjust = 0.5, gp = gpar(col="#091926", fontsize=12)))

#### Figure B2.3: Combined difference plots AMCE vs MM for stealth ####
diff_amce_stealth_rating <- cj(sub_stealth, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~stealthsplit)
diff_mm_stealth_rating <- cj(sub_stealth, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~stealthsplit)
diff_amce_stealth_choice <- cj(sub_stealth, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~stealthsplit)
diff_mm_stealth_choice <- cj(sub_stealth, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"),  by = ~stealthsplit)
diff_amce_stealth_rating$Estimate <- "AMCE (Rating)"
diff_mm_stealth_rating$Estimate <- "Marginal Means (Rating)"
diff_amce_stealth_choice$Estimate <- "AMCE (Choice)"
diff_mm_stealth_choice$Estimate <- "Marginal Means (Choice)"
diff_stealth_rating <- rbind(diff_amce_stealth_rating, diff_mm_stealth_rating)
diff_stealth_choice <- rbind(diff_amce_stealth_choice, diff_mm_stealth_choice)

stealth1 <- plot(diff_stealth_rating, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

stealth2 <- plot(diff_stealth_choice, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

combined_stealth <- ggarrange(stealth1 + theme(axis.ticks.y = element_blank(),
                                               plot.margin = margin(r = 1) ),
                              stealth2 + theme(axis.text.y = element_blank(),
                                               axis.ticks.y = element_blank(),
                                               axis.title.y = element_blank(),
                                               plot.margin = margin(l = 1)),
                              nrow = 1,
                              bottom = textGrob("Estimated differences for 'high stealth attitudes' compared to 'low stealth attitudes'", vjust = 0.5, gp = gpar(col="#091926", fontsize=12)))

#### Figure B2.4: Combined difference plots AMCE vs MM for populism ####
diff_amce_pop_rating <- cj(sub_pop, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~popsplit)
diff_mm_pop_rating <- cj(sub_pop, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~popsplit)
diff_amce_pop_choice <- cj(sub_pop, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~popsplit)
diff_mm_pop_choice <- cj(sub_pop, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability",  "authorization", "format", "consensus", "composition",  "size", "recruitment", "initiative", "issue"), by = ~popsplit)
diff_amce_pop_rating$Estimate <- "AMCE (Rating)"
diff_mm_pop_rating$Estimate <- "Marginal Means (Rating)"
diff_amce_pop_choice$Estimate <- "AMCE (Choice)"
diff_mm_pop_choice$Estimate <- "Marginal Means (Choice)"
diff_pop_rating <- rbind(diff_amce_pop_rating, diff_mm_pop_rating)
diff_pop_choice <- rbind(diff_amce_pop_choice, diff_mm_pop_choice)

pop1 <- plot(diff_pop_rating, feature_headers = FALSE, xlim = c(-0.4,0.5)) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

pop2 <- plot(diff_pop_choice, feature_headers = FALSE, xlim = c(-0.15,0.15)) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

combined_pop <- ggarrange(pop1 + theme(axis.ticks.y = element_blank(),
                                       plot.margin = margin(r = 1) ),
                          pop2 + theme(axis.text.y = element_blank(),
                                       axis.ticks.y = element_blank(),
                                       axis.title.y = element_blank(),
                                       plot.margin = margin(l = 1)),
                          nrow = 1,
                          bottom = textGrob("Estimated differences for 'high populist' compared to 'low populist'", vjust = 0.5, gp = gpar(col="#091926", fontsize=12)))

##### 2. Conditional AMCE choice #####
#### Figure B3.1: Conditional AMCE (facetting) political dissatisfaction ####
mm_by_sat_choice <- cj(sub_satisfaction, Mchoice, id = ~ID, estimate="mm", h0=0.5, feature_labels = feature_names, by = ~satsplit)
mm_by_sat_rating <- cj(sub_satisfaction, Mrating, id = ~ID, estimate="mm", feature_labels = feature_names, by = ~satsplit)
amce_by_sat_choice <- cj(sub_satisfaction, Mchoice, id = ~ID,estimate = "amce",feature_labels = feature_names, by = ~satsplit)
amce_by_sat_rating <- cj(sub_satisfaction, Mrating, id = ~ID,estimate = "amce",feature_labels = feature_names, by = ~satsplit)

diff_amce_sat_cond_choice <- cj(sub_satisfaction, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~satsplit)
amce_sat_choice_diff <- rbind(amce_by_sat_choice, diff_amce_sat_cond_choice)
diff_amce_sat_cond_rating <- cj(sub_satisfaction, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~satsplit)
amce_sat_rating_diff <- rbind(amce_by_sat_rating, diff_amce_sat_cond_rating)

plot(amce_sat_choice_diff, feature_headers = FALSE) +
  ggplot2::facet_grid(feature~BY, scales = "free_y", space="free", switch = "y") +
  theme(strip.background.x = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        strip.text.y = element_blank(),
        panel.spacing.y = unit(1, "lines"),
        panel.border = element_rect(linetype = "blank"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 12),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#304965", "#833d3d"))

#### Figure B3.2: Conditional AMCE (facetting) external efficacy ####
mm_by_external_choice <- cj(sub_external, Mchoice, id = ~ID, estimate = "mm", feature_labels = feature_names, by = ~externalsplit)
mm_by_external_rating <- cj(sub_external, Mrating, id = ~ID, estimate = "mm", feature_labels = feature_names, by = ~externalsplit)
amce_by_external_choice <- cj(sub_external, Mchoice, id = ~ID, estimate="amce", feature_labels = feature_names, by = ~externalsplit)
amce_by_external_rating <- cj(sub_external, Mrating, id = ~ID, estimate="amce", feature_labels = feature_names, by = ~externalsplit)

diff_amce_external_cond_choice <- cj(sub_external, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~externalsplit)
amce_external_choice_diff <- rbind(amce_by_external_choice, diff_amce_external_cond_choice)
diff_amce_external_cond_rating <- cj(sub_external, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~externalsplit)
amce_external_rating_diff <- rbind(amce_by_external_rating, diff_amce_external_cond_rating)

plot(amce_external_choice_diff, feature_headers = FALSE) +
  ggplot2::facet_grid(feature~BY, scales = "free_y", space="free", switch = "y") +
  theme(strip.background.x = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        strip.text.y = element_blank(),
        panel.spacing.y = unit(1, "lines"),
        panel.border = element_rect(linetype = "blank"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 12),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#304965", "#833d3d"))

#### Figure B3.3: Conditional AMCE (facetting) stealth ####
mm_by_stealth_choice <- cj(sub_stealth, Mchoice, id = ~ID, estimate = "mm", feature_labels = feature_names, by = ~stealthsplit)
mm_by_stealth_rating <- cj(sub_stealth, Mrating, id = ~ID, estimate = "mm", feature_labels = feature_names, by = ~stealthsplit)
amce_by_stealth_choice <- cj(sub_stealth, Mchoice, id = ~ID, estimate = "amce", feature_labels = feature_names, by = ~stealthsplit)
amce_by_stealth_rating <- cj(sub_stealth, Mrating, id = ~ID, estimate = "amce", feature_labels = feature_names, by = ~stealthsplit)

diff_amce_stealth_cond_choice <- cj(sub_stealth, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~stealthsplit)
amce_stealth_choice_diff <- rbind(amce_by_stealth_choice, diff_amce_stealth_cond_choice)
diff_amce_stealth_cond_rating <- cj(sub_stealth, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~stealthsplit)
amce_stealth_rating_diff <- rbind(amce_by_stealth_rating, diff_amce_stealth_cond_rating)

plot(amce_stealth_choice_diff, feature_headers = FALSE) +
  ggplot2::facet_grid(feature~BY, scales = "free_y", space="free", switch = "y") +
  theme(strip.background.x = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        strip.text.y = element_blank(),
        panel.spacing.y = unit(1, "lines"),
        panel.border = element_rect(linetype = "blank"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 12),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#304965", "#833d3d"))

#### Figure B3.4: Conditional AMCE (facetting) populism ####
mm_by_pop_choice <- cj(sub_pop, Mchoice, id = ~ID, estimate="mm", h0=0.5, feature_labels = feature_names, by = ~popsplit)
mm_by_pop_rating <- cj(sub_pop, Mrating, id = ~ID, estimate="mm", feature_labels = feature_names, by = ~popsplit)
amce_by_pop_choice <- cj(sub_pop, Mchoice, id = ~ID,estimate = "amce", feature_labels = feature_names, by = ~popsplit)
amce_by_pop_rating <- cj(sub_pop, Mrating, id = ~ID,estimate = "amce", feature_labels = feature_names, by = ~popsplit)

diff_amce_pop_cond_choice <- cj(sub_pop, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~popsplit)
amce_pop_choice_diff <- rbind(amce_by_pop_choice, diff_amce_pop_cond_choice)
diff_amce_pop_cond_rating <- cj(sub_pop, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~popsplit)
amce_pop_rating_diff <- rbind(amce_by_pop_rating, diff_amce_pop_cond_rating)

plot(amce_pop_choice_diff, feature_headers = FALSE) +
  ggplot2::facet_grid(feature~BY, scales = "free_y", space="free", switch = "y") +
  theme(strip.background.x = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        strip.text.y = element_blank(),
        panel.spacing.y = unit(1, "lines"),
        panel.border = element_rect(linetype = "blank"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 12),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#304965", "#833d3d"))

######### Timing variables ####################################################
##### 1. median split #####
mm_by_time_choice <- cj(cj_data, Mchoice, id = ~ID, estimate="mm", h0=0.5, feature_labels = feature_names, by = ~timing)
mm_by_time_rating <- cj(cj_data, Mrating, id = ~ID, estimate="mm", feature_labels = feature_names, by = ~timing)
amce_by_time_choice <- cj(cj_data, Mchoice, id = ~ID,  estimate="amce", feature_labels = feature_names, by = ~timing)
amce_by_time_rating <- cj(cj_data, Mrating, id = ~ID, estimate="amce", feature_labels = feature_names, by = ~timing)

diff_amce_time_choice <- cj(cj_data, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~timing)
amce_time_choice_diff <- rbind(amce_by_time_choice, diff_amce_time_choice)
diff_mm_time_choice <- cj(cj_data, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names, by = ~timing)
mm_time_choice_diff <- rbind(mm_by_time_choice, diff_mm_time_choice)
diff_amce_time_rating <- cj(cj_data, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, by = ~timing)
amce_time_rating_diff <- rbind(amce_by_time_rating, diff_amce_time_rating)
diff_mm_time_rating <- cj(cj_data, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, by = ~timing)
mm_time_rating_diff <- rbind(mm_by_time_rating, diff_mm_time_rating)

#### Figure B4.1: Conditional AMCE (facetting) choice #####
plot(amce_time_choice_diff, feature_headers = FALSE) +
  ggplot2::facet_grid(feature~BY, scales = "free_y", space="free", switch = "y") +
  theme(strip.background.x = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        strip.text.y = element_blank(),
        panel.spacing.y = unit(1, "lines"),
        panel.border = element_rect(linetype = "blank"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 12),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d", "#304965"))

#### Figure B4.2: Conditional AMCE (facetting) median split rating #####
plot(amce_time_rating_diff, feature_headers = FALSE) +
  ggplot2::facet_grid(feature~BY, scales = "free_y", space="free", switch = "y") +
  theme(strip.background.x = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        strip.text.y = element_blank(),
        panel.spacing.y = unit(1, "lines"),
        panel.border = element_rect(linetype = "blank"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 12),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d", "#304965"))

##### 2. Quantiles #####
quantile(cj_data$page_p_screen1_timing, probs = c(.25, .5, .75))
quantile(cj_data$page_p_screen2_timing, probs = c(.25, .5, .75))
quantile(cj_data$page_p_screen3_timing, probs = c(.25, .5, .75))
quantile(cj_data$page_p_screen4_timing, probs = c(.25, .5, .75))
quantile(cj_data$page_p_screen5_timing, probs = c(.25, .5, .75))
quantile(cj_data$page_p_screen6_timing, probs = c(.25, .5, .75))

# 1st quantile
cj_data$select_timing1a <- ifelse(cj_data$comparision_table == 1 & cj_data$page_p_screen1_timing > 31.08, 1, 0)
cj_data$select_timing2a <- ifelse(cj_data$comparision_table == 2 & cj_data$page_p_screen2_timing > 19.87, 1, 0)
cj_data$select_timing3a <- ifelse(cj_data$comparision_table == 3 & cj_data$page_p_screen3_timing > 17.45, 1, 0)
cj_data$select_timing4a <- ifelse(cj_data$comparision_table == 4 & cj_data$page_p_screen4_timing > 16.29, 1, 0)
cj_data$select_timing5a <- ifelse(cj_data$comparision_table == 5 & cj_data$page_p_screen5_timing > 14.99, 1, 0)
cj_data$select_timing6a <- ifelse(cj_data$comparision_table == 6 & cj_data$page_p_screen6_timing > 14.35, 1, 0)
cj_data$timing_a <- rowSums(subset(cj_data, select = c(select_timing1a, select_timing2a, select_timing3a, select_timing4a, select_timing5a, select_timing6a)), na.rm = TRUE)
cj_data$timing_a <- factor(cj_data$timing_a, labels = c("below 1st quantile", "above 1st quantile"))

#3rd quantile
cj_data$select_timing1b <- ifelse(cj_data$comparision_table == 1 & cj_data$page_p_screen1_timing > 74.13, 1, 0)
cj_data$select_timing2b <- ifelse(cj_data$comparision_table == 2 & cj_data$page_p_screen2_timing > 49.09, 1, 0)
cj_data$select_timing3b <- ifelse(cj_data$comparision_table == 3 & cj_data$page_p_screen3_timing > 43.87, 1, 0)
cj_data$select_timing4b <- ifelse(cj_data$comparision_table == 4 & cj_data$page_p_screen4_timing > 39.15, 1, 0)
cj_data$select_timing5b <- ifelse(cj_data$comparision_table == 5 & cj_data$page_p_screen5_timing > 37.02, 1, 0)
cj_data$select_timing6b <- ifelse(cj_data$comparision_table == 6 & cj_data$page_p_screen6_timing > 35.23, 1, 0)
cj_data$timing_b <- rowSums(subset(cj_data, select = c(select_timing1b, select_timing2b, select_timing3b, select_timing4b, select_timing5b, select_timing6b)), na.rm = TRUE)
cj_data$timing_b <- factor(cj_data$timing_b, labels = c("below 3rd quantile", "above 3rd quantile"))

#### Figure B4.3: Differences for engagement with scenarios AMCE choice ####
diff_amce_1st_choice <- cj(cj_data, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing_a)
diff_amce_2nd_choice <- cj(cj_data, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing)
diff_amce_3rd_choice <- cj(cj_data, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing_b)

diff_engagement_choice <-bind_rows(
  diff_amce_1st_choice %>% mutate(model = "above 1st quantile"),
  diff_amce_2nd_choice %>% mutate(model = "above median"),
  diff_amce_3rd_choice %>% mutate(model = "above 3rd quantile"))

diff_engagement_choice <- diff_engagement_choice %>% select(BY:level, estimate:model)
diff_engagement_choice$model <- factor(diff_engagement_choice$model, levels = c("above 1st quantile", "above median", "above 3rd quantile"))

plot(diff_engagement_choice, feature_headers = FALSE, xlab = "Estimated Difference (AMCE)") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 10),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = rep("#833d3d", 9))

#### Figure B4.4: Differences for engagement with scenarios MM choice ####
diff_mm_1st_choicemm <- cj(cj_data, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing_a)
diff_mm_2nd_choicemm <- cj(cj_data, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing)
diff_mm_3rd_choicemm <- cj(cj_data, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing_b)

diff_engagement_choicemm <-bind_rows(
  diff_mm_1st_choicemm %>% mutate(model = "above 1st quantile"),
  diff_mm_2nd_choicemm %>% mutate(model = "above median"),
  diff_mm_3rd_choicemm %>% mutate(model = "above 3rd quantile"))

diff_engagement_choicemm <- diff_engagement_choicemm %>% select(BY:upper, model)
diff_engagement_choicemm$model <- factor(diff_engagement_choicemm$model, levels = c("above 1st quantile", "above median", "above 3rd quantile"))

plot(diff_engagement_choicemm, feature_headers = FALSE, xlab = "Estimated Difference (MM)") + 
  ggplot2::facet_wrap(~model, ncol = 4) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 10),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = rep("#833d3d", 9))

#### Figure B4.5: Differences for engagement with scenarios AMCE rating ####
diff_amce_1st_rating <- cj(cj_data, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing_a)
diff_amce_2nd_rating <- cj(cj_data, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing)
diff_amce_3rd_rating <- cj(cj_data, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing_b)

diff_engagement_rating <-bind_rows(
  diff_amce_1st_rating %>% mutate(model = "above 1st quantile"),
  diff_amce_2nd_rating %>% mutate(model = "above median"),
  diff_amce_3rd_rating %>% mutate(model = "above 3rd quantile"))

diff_engagement_rating <- diff_engagement_rating %>% select(BY:level, estimate:model)
diff_engagement_rating$model <- factor(diff_engagement_rating$model, levels = c("above 1st quantile", "above median", "above 3rd quantile"))

plot(diff_engagement_rating, feature_headers = FALSE, xlab = "Estimated Difference (AMCE)") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 10),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = rep("#833d3d", 9))

#### Figure B4.6: Differences for engagement with scenarios MM rating ####
diff_mm_1st_ratingmm <- cj(cj_data, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing_a)
diff_mm_2nd_ratingmm <- cj(cj_data, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing)
diff_mm_3rd_ratingmm <- cj(cj_data, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~timing_b)

diff_engagement_ratingmm <-bind_rows(
  diff_mm_1st_ratingmm %>% mutate(model = "above 1st quantile"),
  diff_mm_2nd_ratingmm %>% mutate(model = "above median"),
  diff_mm_3rd_ratingmm %>% mutate(model = "above 3rd quantile"))

diff_engagement_ratingmm <- diff_engagement_ratingmm %>% select(BY:upper, model)
diff_engagement_ratingmm$model <- factor(diff_engagement_ratingmm$model, levels = c("above 1st quantile", "above median", "above 3rd quantile"))

plot(diff_engagement_ratingmm, feature_headers = FALSE, xlab = "Estimated Difference (MM)") + 
  ggplot2::facet_wrap(~model, ncol = 4) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#091926", size = 10),
        axis.title = element_text(colour = "#091926", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = rep("#833d3d", 9))

######### Experience ##########################################################
##### 1. Prepare subsets #####
cj_experience <- cj_data %>% select(ID:authorization, q21, outcome.favorability) %>% drop_na()

# 1. experienced = more than heard about 
cj_experience$experience <- ifelse(cj_experience$q21 > 2, 1, 0)
cj_experience$experience <- factor(cj_experience$experience, labels = c("inexperienced", "experienced"))
# 2. experienced = including heard about 
cj_experience$experience2 <- ifelse(cj_experience$q21 > 1, 1, 0)
cj_experience$experience2 <- factor(cj_experience$experience2, labels = c("inexperienced", "experienced"))

# positive and negative experiences
sub_experience <- cj_data %>% select(ID:authorization, q22, outcome.favorability) %>% drop_na()
sub_experience$experience[sub_experience$q22 >= 5] <- 1L
sub_experience$experience[sub_experience$q22 < 5] <- 2L
sub_experience$experience <- factor(sub_experience$experience, 1:2, c("positive experiences", "negative experiences"))

##### 2. Combined: Difference plots AMCE vs MM experiences #####
#### Figure 5.1: Difference plots for experiences with DCFs ####
diff_amce_exp1_rating <- cj(cj_experience, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~experience2)
diff_mm_exp1_rating <- cj(cj_experience, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~experience2)
diff_amce_exp1_choice <- cj(cj_experience, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~experience2)
diff_mm_exp1_choice <- cj(cj_experience, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~experience2)
diff_amce_exp1_rating$Estimate <- "AMCE (Rating)"
diff_mm_exp1_rating$Estimate <- "Marginal Means (Rating)"
diff_amce_exp1_choice$Estimate <- "AMCE (Choice)"
diff_mm_exp1_choice$Estimate <- "Marginal Means (Choice)"
diff_exp1_rating <- rbind(diff_amce_exp1_rating, diff_mm_exp1_rating)
diff_exp1_choice <- rbind(diff_amce_exp1_choice, diff_mm_exp1_choice)

expA1 <- plot(diff_exp1_rating, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

expA2 <- plot(diff_exp1_choice, feature_headers = FALSE, xlim = c(-0.1,0.1)) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

combined_exp1 <- ggarrange(expA1 + theme(axis.ticks.y = element_blank(),
                                         plot.margin = margin(r = 1) ),
                           expA2 + theme(axis.text.y = element_blank(),
                                         axis.ticks.y = element_blank(),
                                         axis.title.y = element_blank(),
                                         plot.margin = margin(l = 1)),
                           nrow = 1,
                           bottom = textGrob("Estimated differences for 'at least some experiences' compared to 'no experiences'", vjust = 0.5, gp = gpar(col="#091926", fontsize=12)))

#### Figure 5.2: Difference plots for experiences with DCFs Positive / negative ####
diff_amce_exp2_rating <- cj(sub_experience, Mrating, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~experience)
diff_mm_exp2_rating <- cj(sub_experience, Mrating, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~experience)
diff_amce_exp2_choice <- cj(sub_experience, Mchoice, id = ~ID, estimate="amce_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~experience)
diff_mm_exp2_choice <- cj(sub_experience, Mchoice, id = ~ID, estimate="mm_diff", feature_labels = feature_names, feature_order = c("outcome.favorability", "authorization", "format", "consensus", "composition", "size", "recruitment", "initiative", "issue"), by = ~experience)
diff_amce_exp2_rating$Estimate <- "AMCE (Rating)"
diff_mm_exp2_rating$Estimate <- "Marginal Means (Rating)"
diff_amce_exp2_choice$Estimate <- "AMCE (Choice)"
diff_mm_exp2_choice$Estimate <- "Marginal Means (Choice)"
diff_exp2_rating <- rbind(diff_amce_exp2_rating, diff_mm_exp2_rating)
diff_exp2_choice <- rbind(diff_amce_exp2_choice, diff_mm_exp2_choice)


expB1 <- plot(diff_exp2_rating, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

expB2 <- plot(diff_exp2_choice, feature_headers = FALSE) + 
  ggplot2::facet_wrap(~Estimate, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "#515151", size = 12),
        axis.text.x = element_text(colour = "#515151", size = 12),
        axis.title = element_blank(),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "#091926", size = 14),
        legend.position = "none")+
  geom_point(size=2) +
  scale_color_manual(values = c("#833d3d"))

combined_exp2 <- ggarrange(expB1 + theme(axis.ticks.y = element_blank(),
                                         plot.margin = margin(r = 1) ),
                           expB2 + theme(axis.text.y = element_blank(),
                                         axis.ticks.y = element_blank(),
                                         axis.title.y = element_blank(),
                                         plot.margin = margin(l = 1)),
                           nrow = 1,
                           bottom = textGrob("Estimated differences for 'negative experiences' compared to 'positive experiences'", vjust = 0.5, gp = gpar(col="#091926", fontsize=12)))
